# -*- coding: utf-8 -*-
"""
The python script to generate Figure 3. 
It includes 1. read monthly air temperature, huv, sensible heat flux, and 
latent heat flux data on three pairs of land tiles
2. Calculate the Bowen ratio from sensible and latent heat flux
3. Calculate the monthly air temperature, huv, and Bowen ratio in the 
North hemisphere for three pairs of land tiles
3. Plot the seasonal variations of air temperature, huv, and Bowen ratio for
three pairs of land tiles 

"""
import netCDF4 as nc
import numpy as np
import matplotlib.pyplot as plt
import warnings
warnings.filterwarnings("ignore")

Dir ='G:\\SSP5UC_DataPreparation\\deposit\\Monthly\\' # The directory you are working with
LonNum=288 # Number of longitudinal grids
LatNum=192 # Number of latitudinal grids

#Read the area of each grid
AreaFile=nc.Dataset('G:\\SSP5UC_DataPreparation\\deposit\\Surface\\b.e21.SSP585UC_IndividualSoil_surfdata_0.9x1.25_16pfts.nc')
area=AreaFile.variables['AREA'][96:192,:]

def AreaWeightedMean(VarData):
    """
    AreaWeightedMean
    Description: a function to calculate the area-weighed mean value of a variable
    Input： VarData--a 2D data in the shape of (96,288)
    Output: GlobalMean--area-weighed global mean value of VarData
     
    """  
    MaskedArea=np.ma.masked_where(np.ma.getmask(VarData[0]),area)
    GlobalMean=np.nansum(VarData*(MaskedArea/np.sum(MaskedArea)),axis=(1,2))
    return GlobalMean

def ReadVar(Varname,Perturbed,Base):
    """
    ReadVar
    Description: a function to read monthly state varibles for the perturbed and base land tiles
                 and then calculate their monthly North hemisphere(NH) mean values
    Input： Varname--the short for the targeted variable (tas/huv)
    Perturbed--the name of the perturned land tiles (u_Urban/l_Lake/g_Grass)
    Base--the name of the base land tiles (r_Rural/r_Rural/t_Tree)
    Output: VarpMean--monthly mean values in the North hemisphere for the perturbed land tile
            VarbMean--monthly mean values in the North hemisphere for the base land tile
     
    """  
    FileP=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Perturbed+'_Monthly_201501-210012.nc')
    FileB=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Base+'_Monthly_201501-210012.nc')
    Varp=FileP.variables[Varname+Perturbed[0]][0:360,96:192,:] # Here 96:192 means only read data for NH
    Varp=np.nanmean(np.ma.array(np.split(Varp, len(Varp[:,0,0])/12)),axis=0) #calculate mean monthly result during 2015-2044
    Varb=FileB.variables[Varname+Base[0]][0:360,96:192,:]
    Varb=np.nanmean(np.ma.array(np.split(Varb, len(Varb[:,0,0])/12)),axis=0)
    Varp=np.ma.masked_where(np.ma.getmask(MaskArray),Varp) # Only keep grids with both perturbed and base land tiles
    Varb=np.ma.masked_where(np.ma.getmask(MaskArray),Varb) # Only keep grids with both perturbed and base land tiles
    VarpMean=AreaWeightedMean(Varp);VarbMean=AreaWeightedMean(Varb);
    return VarpMean,VarbMean

def ReadVarNoMean(Varname,Perturbed,Base):
    """
    ReadVar
    Description: a function to read monthly state varibles for the perturbed and base land tiles

    Input： Varname--the short for the targeted variable (hfls/hfss)
    Perturbed--the name of the perturned land tiles (u_Urban/l_Lake/g_Grass)
    Base--the name of the base land tiles (r_Rural/r_Rural/t_Tree)
    Output: Varp--monthly data of the targeted variable in NH during 2015-2044 for the perturbed land tile
            Varb--monthly data of the targeted variable in NH during 2015-2044 for the base land tile
    """  
    FileP=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Perturbed+'_Monthly_201501-210012.nc')
    FileB=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_'+Varname+Base+'_Monthly_201501-210012.nc')
    Varp=FileP.variables[Varname+Perturbed[0]][0:360,96:192,:] # Here 96:192 means only read data for NH
    Varp=np.nanmean(np.ma.array(np.split(Varp, len(Varp[:,0,0])/12)),axis=0) #calculate mean monthly result during 2015-2044
    Varb=FileB.variables[Varname+Base[0]][0:360,96:192,:]
    Varb=np.nanmean(np.ma.array(np.split(Varb, len(Varb[:,0,0])/12)),axis=0)
    Varp=np.ma.masked_where(np.ma.getmask(MaskArray),Varp) # Only keep grids with both perturbed and base land tiles
    Varb=np.ma.masked_where(np.ma.getmask(MaskArray),Varb) # Only keep grids with both perturbed and base land tiles
    return Varp,Varb

def MainRead(Perturbedname,Basename):
    """
    MainRead
    Description: a main function to calculate monthly mean values of tas, huv, and Bowen ratio 
                in the North hemisphere for the perturbed and base land tile by calling ReadVar and ReadVarNoMean 

    Input： Perturbedname--the name of the perturned land tiles (u_Urban/l_Lake/g_Grass)
           Basename--the name of the base land tiles (r_Rural/r_Rural/t_Tree)
    Output: tasp--monthly mean air temperature in the North hemisphere for the perturbed land tile
    tasb--monthly mean air temperature in the North hemisphere for the base land tile
    huvp--monthly mean vapor pressure in the North hemisphere for the perturbed land tile
    huvb--monthly mean vapor pressure in the North hemisphere for the base land tile
    PerturbedBowenRatio--monthly mean Bowen ratio in the North hemisphere for the perturbed land tile
    BaseBowenRatio--monthly mean Bowen ratio in the North hemisphere for the base land tile
    """ 
    [tasp,tasb]=ReadVar('tas',Perturbedname,Basename)
    [huvp,huvb]=ReadVar('huv',Perturbedname,Basename)
    [hflsp,hflsb]=ReadVar('hfls',Perturbedname,Basename)
    [hfssp,hfssb]=ReadVar('hfss',Perturbedname,Basename)
    [hflspNoMean,hflsbNoMean]=ReadVarNoMean('hfls',Perturbedname,Basename)
    [hfsspNoMean,hfssbNoMean]=ReadVarNoMean('hfss',Perturbedname,Basename)
    #BowenRatio
    BowenRatioPerturbed = hfsspNoMean/hflspNoMean
    BowenRatioBase = hfssbNoMean/hflsbNoMean
    # When the sensible and latent heat flux are too small (<0.15 W/m2) or the Bowen ratio is too larger(>25), the result is not quite meaningful
    PerturbedBowenRatio = np.ma.masked_where((np.absolute(hflspNoMean)<0.15)|(np.absolute(hfsspNoMean)<0.15)|(np.absolute(BowenRatioPerturbed)>25),BowenRatioPerturbed)
    BaseBowenRatio =np.ma.masked_where((np.absolute(hflsbNoMean)<0.15)|(np.absolute(hfssbNoMean)<0.15)|(np.absolute(BowenRatioBase)>25),BowenRatioBase) #
    # Make sure the unmasked grids are consistent for the perturbed and base land tile
    PerturbedBowenRatio = np.ma.masked_where(np.ma.getmask(BaseBowenRatio),PerturbedBowenRatio)
    BaseBowenRatio= np.ma.masked_where(np.ma.getmask(PerturbedBowenRatio),BaseBowenRatio)
    
    PerturbedBowenRatio=AreaWeightedMean(PerturbedBowenRatio)
    BaseBowenRatio=AreaWeightedMean(BaseBowenRatio)
    return tasp,tasb,huvp,huvb,PerturbedBowenRatio,BaseBowenRatio


#--------------------------Urban vs Rural-----------------------------
PerturbednameInput='u_Urban';BasenameInput='r_Rural';
# Read two arrays from the perturbed and base land tiles
# Their maskes are useful to ensure only grids with both perturbed and base land tiles are used
MaskU=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+PerturbednameInput+'_Monthly_201501-210012.nc')
MaskR=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+BasenameInput+'_Monthly_201501-210012.nc')
MaskArray=MaskU.variables['tas'+PerturbednameInput[0]][0:12,96:192,:]-MaskR.variables['tas'+BasenameInput[0]][0:12,96:192,:]
[tasUrban,tasRural,huvUrban,huvRural,BRUrban,BRRural]=MainRead(PerturbednameInput,BasenameInput)

#--------------------------Lake vs Rural-----------------------------
PerturbednameInput='l_Lake';BasenameInput='r_Rural';
# Read two arrays from the perturbed and base land tiles
# Their maskes are useful to ensure only grids with both perturbed and base land tiles are used
MaskU=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+PerturbednameInput+'_Monthly_201501-210012.nc')
MaskR=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+BasenameInput+'_Monthly_201501-210012.nc')
MaskArray=MaskU.variables['tas'+PerturbednameInput[0]][0:12,96:192,:]-MaskR.variables['tas'+BasenameInput[0]][0:12,96:192,:]
[tasLake,tasNonlake,huvLake,huvNonlake,BRLake,BRNonlake]=MainRead(PerturbednameInput,BasenameInput)

#--------------------------Grass vs Tree-----------------------------
PerturbednameInput='g_Grass';BasenameInput='t_Tree'#
# Read two arrays from the perturbed and base land tiles
# Their maskes are useful to ensure only grids with both perturbed and base land tiles are used
MaskU=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+PerturbednameInput+'_Monthly_201501-210012.nc')
MaskR=nc.Dataset(Dir+'b.e21.SSP585UC_IndividualSoil_tas'+BasenameInput+'_Monthly_201501-210012.nc')
MaskArray=MaskU.variables['tas'+PerturbednameInput[0]][0:12,96:192,:]-MaskR.variables['tas'+BasenameInput[0]][0:12,96:192,:]
[tasGrass,tasTree,huvGrass,huvTree,BRGrass,BRTree]=MainRead(PerturbednameInput,BasenameInput)

# Define colors for three land tile pairs
UrbanColor='#CC0000';RuralColor='#16A085'
GrassColor='sandybrown';TreeColor='forestgreen'
LakeColor='cornflowerblue';NonLakeColor='salmon'

fig,ax = plt.subplots(nrows=3, ncols=3, figsize=(9,4.7))
Months=np.arange(1,13,1)

ax[0,0].plot(Months,tasUrban[:]-273.15,color=UrbanColor,linewidth=1.4,label="Urban Tas")
ax[0,0].plot(Months,tasRural[:]-273.15,color=RuralColor,linewidth=1.4,label="Rural Tas")
# ax[0,0].set_xlabel('Month',fontsize=14,labelpad=1)
# plt.ylabel('$\it{Ta}$'+' (°C)',fontsize=11,labelpad=1)
# ax[0,0].set_ylabel('tas (°C)',fontsize=11,labelpad=1)
ax[0,0].set_xticks(np.arange(0, 14, 2))
ax[0,0].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[0,0].tick_params(axis="both", labelsize=11,direction='in')
ax[0,0].set_xlim((-0.5,12.5))
ax[0,0].set_ylim((4,28))
ax[0,0].set_yticks(np.arange(8, 27, 6))
ax[0,0].text(0, 24, '(a)', horizontalalignment='left',verticalalignment='bottom', size=12)

ax0 = ax[0,0].twinx()
ax0.plot(Months, tasUrban[:]-tasRural[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Urban-Rural Ta")
ax0.tick_params(axis="both", labelsize=10,direction='in')
ax0.set_ylim((0.5,1))
ax0.set_yticks(np.arange(0.6, 0.91, 0.1))

# ax[0,1].plot(Months,GridMeanVAPOR_PRES[:],'k',linewidth=1.4,label="Grid vapore pressure")
ax[0,1].plot(Months,huvUrban[:]*0.01,color=UrbanColor,linewidth=1.4,label="Urban huv")
ax[0,1].plot(Months,huvRural[:]*0.01,color=RuralColor,linewidth=1.4,label="Rural huv")
# ax[0,1].set_xlabel('Month',fontsize=14,labelpad=1)
# ax[0,1].set_ylabel('$\it{Ea}$'+' (hPa)',fontsize=11,labelpad=1)
# ax[0,1].set_ylabel('huv (hPa)',fontsize=11,labelpad=1)
ax[0,1].set_xticks(np.arange(0, 14, 2))
ax[0,1].set_yticks(np.arange(9, 23, 6))
ax[0,1].set_ylim( (5,25) )
ax[0,1].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[0,1].tick_params(axis="both", labelsize=11,direction='in')
ax[0,1].set_xlim((-0.5,12.5))
ax[0,1].text(0, 21.5, '(b)', horizontalalignment='left',verticalalignment='bottom', size=12)
# ax[0,1].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax1 = ax[0,1].twinx()
ax1.plot(Months, (huvUrban[:]-huvRural[:])*0.01,color='dimgrey',linestyle='dashed',linewidth=0.8,label="Urban-Rural huv")
ax1.tick_params(axis="both", labelsize=10,direction='in')
ax1.set_ylim((-0.7,0))
ax1.set_yticks(np.arange(-0.6, -0.1, 0.2))

# ax[0,1].plot(Months,GridBowenRatio[:],'k',linewidth=1.4,label="Grid Bowen Ratio")
ax[0,2].plot(Months,BRUrban[:],color=UrbanColor,linewidth=1.4,label="Urban")
ax[0,2].plot(Months,BRRural[:],color=RuralColor,linewidth=1.4,label="Rural")
# ax[0,2].set_xlabel('Month',fontsize=14,labelpad=1)
# ax[0,2].set_ylabel('Bowen Ratio',fontsize=11,labelpad=1)
ax[0,2].set_xticks(np.arange(0, 14, 2))
ax[0,2].set_ylim( (0, 3.5) )
ax[0,2].set_yticks(np.arange(0.5, 3, 1))

ax[0,2].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[0,2].tick_params(axis="both", labelsize=11,direction='in')
ax[0,2].set_xlim((-0.5,12.5))
ax[0,2].text(0, 2.95, '(c)', horizontalalignment='left',verticalalignment='bottom', size=12)
#ax[0,1].ylim( (7.5, 8.5) )
ax2 = ax[0,2].twinx()
ax2.plot(Months, BRUrban[:]-BRRural[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Urban - Rural")
ax2.set_yticks(np.arange(1.1, 2.4, 0.5))
ax2.tick_params(axis="both", labelsize=10,direction='in')
lines, labels = ax[0,2].get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2,bbox_to_anchor=(1.23,0.24),frameon=False, ncol=1,loc='lower left', prop={'size': 7})
#0.56
# ax[0,2].legend(bbox_to_anchor=(-1.8,-2.8),frameon=True, ncol=2,loc='lower left', prop={'size': 11}) #, ncol=3
# ax2.legend(bbox_to_anchor=(-1.8,-2.8),frameon=True, ncol=1,loc='lower left', prop={'size': 11})

ax[1,0].plot(Months,tasGrass[:]-273.15,color=GrassColor,linewidth=1.4,label="Grass Tas")
ax[1,0].plot(Months,tasTree[:]-273.15,color=TreeColor,linewidth=1.4,label="Tree Tas")
# ax[1,0].set_xlabel('Month',fontsize=14,labelpad=1)
# plt.ylabel('$\it{Ta}$'+' (°C)',fontsize=11,labelpad=1)
ax[1,0].set_ylabel('tas (°C)',fontsize=13,labelpad=3,fontweight='bold')
ax[1,0].set_xticks(np.arange(0, 14, 2))
ax[1,0].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[1,0].tick_params(axis="both", labelsize=11,direction='in')
ax[1,0].set_xlim((-0.5,12.5))
ax[1,0].set_ylim( (0, 25) )
ax[1,0].set_yticks(np.arange(5, 21, 5))
# ax[1,0].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=3,loc='lower left', prop={'size': 6}) #
ax[1,0].text(0, 21.2, '(d)', horizontalalignment='left',verticalalignment='bottom', size=12)

ax0 = ax[1,0].twinx()
ax0.plot(Months, tasGrass[:]-tasTree[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Grass-Tree Ta")
ax0.tick_params(axis="both", labelsize=10,direction='in')

# ax[1,1].plot(Months,GridMeanVAPOR_PRES[:],'k',linewidth=1.4,label="Grid vapore pressure")
ax[1,1].plot(Months,huvGrass[:]*0.01,color=GrassColor,linewidth=1.4,label="Grass huv")
ax[1,1].plot(Months,huvTree[:]*0.01,color=TreeColor,linewidth=1.4,label="Tree huv")
# ax[1,1].set_xlabel('Month',fontsize=14,labelpad=1)
# ax[1,1].set_ylabel('$\it{Ea}$'+' (hPa)',fontsize=11,labelpad=1)
ax[1,1].set_ylabel('huv (hPa)',fontsize=13,labelpad=3,fontweight='bold')
ax[1,1].set_xticks(np.arange(0, 14, 2))
ax[1,1].set_yticks(np.arange(6, 21, 4))
ax[1,1].set_ylim( (3, 20) )
ax[1,1].set_xlim((-0.5,12.5))

ax[1,1].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[1,1].tick_params(axis="both", labelsize=11,direction='in')
ax[1,1].text(0, 17, '(e)', horizontalalignment='left',verticalalignment='bottom', size=12)

ax1 = ax[1,1].twinx()
ax1.plot(Months, (huvGrass[:]-huvTree[:])*0.01,color='dimgrey',linestyle='dashed',linewidth=0.8,label="Grass-Tree huv")
ax1.set_ylim( (-0.03, 0.07) )
ax1.set_yticks(np.arange(-0.01, 0.05, 0.03))
ax1.tick_params(axis="both", labelsize=10,direction='in')

# ax[1,1].plot(Months,GridBowenRatio[:],'k',linewidth=1.4,label="Grid Bowen Ratio")
ax[1,2].plot(Months,BRGrass[:],color=GrassColor,linewidth=1.4,label="Grass")
ax[1,2].plot(Months,BRTree[:],color=TreeColor,linewidth=1.4,label="Tree")
# ax[1,2].set_xlabel('Month',fontsize=14,labelpad=1)
ax[1,2].set_ylabel('Bowen ratio',fontsize=12,labelpad=3,fontweight='bold')
ax[1,2].set_xticks(np.arange(0, 14, 2))
ax[1,2].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[1,2].tick_params(axis="both", labelsize=11,direction='in')
ax[1,2].set_xlim((-0.5,12.5))
ax[1,2].set_ylim( (0, 3.9) )
ax[1,2].set_yticks(np.arange(1, 3.1, 1))
ax[1,2].text(0, 3.23, '(f)', horizontalalignment='left',verticalalignment='bottom', size=12)

#ax[1,1].ylim( (7.5, 8.5) )
# ax[1,2].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax2 = ax[1,2].twinx()
ax2.plot(Months, BRGrass[:]-BRTree[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Grass - Tree")
ax2.tick_params(axis="both", labelsize=10,direction='in')
ax2.set_ylim( (-1.7, 0.7) )
ax2.set_yticks(np.arange(-1.5, 0.4, 0.8))

lines, labels = ax[1,2].get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2,bbox_to_anchor=(1.23,0.24),frameon=False, ncol=1,loc='lower left', prop={'size': 7})


ax[2,0].plot(Months,tasLake[:]-273.15,color=LakeColor,linewidth=1.4,label="Lake Tas")
ax[2,0].plot(Months,tasNonlake[:]-273.15,color=NonLakeColor,linewidth=1.4,label="Nonlake Tas")
ax[2,0].set_xlabel('Month',fontsize=14,labelpad=1)
# plt.ylabel('$\it{Ta}$'+' (°C)',fontsize=11,labelpad=1)
# ax[2,0].set_ylabel('tas (°C)',fontsize=11,labelpad=1)
ax[2,0].set_xticks(np.arange(0, 14, 2))
# ax[2,0].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[2,0].tick_params(axis="both", labelsize=11,direction='in')
ax[2,0].set_xlim((-0.5,12.5))
ax[2,0].set_yticks(np.arange(0, 21, 10))
ax[2,0].set_ylim( (-7, 25) )
ax[2,0].text(0, 18.7, '(g)', horizontalalignment='left',verticalalignment='bottom', size=12)

#plt.ylim( (-5, 25) )
# ax[2,0].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=3,loc='lower left', prop={'size': 6}) #

ax0 = ax[2,0].twinx()
ax0.plot(Months, tasLake[:]-tasNonlake[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Lake-Nonlake Ta")
ax0.set_yticks(np.arange(-0.2, 0.3, 0.2))
ax0.set_ylim( (-0.3, 0.5) )
ax0.tick_params(axis="both", labelsize=10,direction='in')

# ax[2,1].plot(Months,GridMeanVAPOR_PRES[:],'k',linewidth=1.4,label="Grid vapore pressure")
ax[2,1].plot(Months,huvLake[:]*0.01,color=LakeColor,linewidth=1.4,label="Lake huv")
ax[2,1].plot(Months,huvNonlake[:]*0.01,color=NonLakeColor,linewidth=1.4,label="Nonlake huv")
ax[2,1].set_xlabel('Month',fontsize=14,labelpad=1)
# ax[2,1].set_ylabel('$\it{Ea}$'+' (hPa)',fontsize=11,labelpad=1)
# ax[2,1].set_ylabel('huv (hPa)',fontsize=11,labelpad=1)
ax[2,1].set_xticks(np.arange(0, 14, 2))
ax[2,1].set_yticks(np.arange(5,19, 6))
ax[2,1].set_ylim( (0, 20) )
ax[2,1].set_xlim((-0.5,12.5))
ax[2,1].text(0, 16.5, '(h)', horizontalalignment='left',verticalalignment='bottom', size=12)

# ax[2,1].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[2,1].tick_params(axis="both", labelsize=11,direction='in')
# ax[2,1].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax1 = ax[2,1].twinx()
ax1.plot(Months, (huvLake[:]-huvNonlake[:])*0.01,color='dimgrey',linestyle='dashed',linewidth=0.8,label="Lake-Nonlake huv")
ax1.tick_params(axis="both", labelsize=10,direction='in')
ax1.set_yticks(np.arange(0.4,1.1, 0.3))
ax1.set_ylim( (0.2, 1.2) )

# ax[2,1].plot(Months,GridBowenRatio[:],'k',linewidth=1.4,label="Grid Bowen Ratio")
ax[2,2].plot(Months,BRLake[:],color=LakeColor,linewidth=1.4,label="Lake")
ax[2,2].plot(Months,BRNonlake[:],color=NonLakeColor,linewidth=1.4,label="Nonlake")
ax[2,2].set_xlabel('Month',fontsize=14,labelpad=1)
# ax[2,2].set_ylabel('Bowen Ratio',fontsize=11,labelpad=1)
ax[2,2].set_xticks(np.arange(0, 14, 2))
# ax[2,2].set_xticklabels([ ]*len(np.arange(0, 14, 2)),fontsize=12)
ax[2,2].tick_params(axis="both", labelsize=11,direction='in')
ax[2,2].set_xlim((-0.5,12.5))
ax[2,2].set_ylim( (-0.5, 3.2) )
ax[2,2].set_yticks(np.arange(0, 2.2, 1))
ax[2,2].text(0, 2.6, '(i)', horizontalalignment='left',verticalalignment='bottom', size=12)
#ax[2,1].ylim( (7.5, 8.5) )
# ax[2,2].legend(bbox_to_anchor=(0,0.98),frameon=False, ncol=2,loc='lower left', prop={'size': 5.7}) #, ncol=3
ax2 = ax[2,2].twinx()
ax2.plot(Months, BRLake[:]-BRNonlake[:],color='dimgrey',linestyle='dashed',linewidth=0.8,label="Lake - Nonlake")
ax2.tick_params(axis="both", labelsize=10,direction='in')
ax2.set_yticks(np.arange(-2.2,-0.5, 0.8))
ax2.set_ylim( (-3, 0) )

lines, labels = ax[2,2].get_legend_handles_labels()
lines2, labels2 = ax2.get_legend_handles_labels()
ax2.legend(lines + lines2, labels + labels2,bbox_to_anchor=(1.23,0.24),frameon=False, ncol=1,loc='lower left', prop={'size': 7})

plt.subplots_adjust(top=0.95,
bottom=0.19,
left=0.096,
right=0.805,
hspace=0.055,
wspace=0.605) 
# plt.savefig('E:\\Keer_work\\DataPaper\\Figure3.png', dpi=450)
